{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian Data Analysis, 3rd ed\n",
    "##  Chapter 3, demo 5\n",
    "\n",
    "Demonstrate a normal model for the Newcomb's data (BDA3 p. 66)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os, sys\n",
    "# add utilities directory to path\n",
    "util_path = os.path.abspath(os.path.join(os.path.pardir, 'utilities_and_data'))\n",
    "if util_path not in sys.path and os.path.exists(util_path):\n",
    "    sys.path.insert(0, util_path)\n",
    "\n",
    "# import from utilities\n",
    "import plot_tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# edit default plot settings\n",
    "plt.rc('font', size=12)\n",
    "# apply custom background plotting style\n",
    "plt.style.use(plot_tools.custom_styles['gray_background'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# data\n",
    "data_path = os.path.abspath(\n",
    "    os.path.join(\n",
    "        os.path.pardir,\n",
    "        'utilities_and_data',\n",
    "        'light.txt'\n",
    "    )\n",
    ")\n",
    "y = np.loadtxt(data_path)\n",
    "# sufficient statistics\n",
    "n = len(y)\n",
    "s2 = np.var(y, ddof=1)  # ddof=1 -> sample estimate\n",
    "my = np.mean(y)\n",
    "\n",
    "# filtered data\n",
    "y_pos = y[y > 0]\n",
    "# sufficient statistics\n",
    "n_pos = len(y_pos)\n",
    "s2_pos = np.var(y_pos, ddof=1)\n",
    "my_pos = np.mean(y_pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# for mu, compute the density in these points\n",
    "tl1 = [10, 40]\n",
    "t1 = np.linspace(tl1[0], tl1[1], 100)\n",
    "\n",
    "# compute the exact marginal density for mu\n",
    "# multiplication by 1./sqrt(s2/n) is due to the transformation of variable\n",
    "# z=(x-mean(y))/sqrt(s2/n), see BDA3 p. 21\n",
    "pm_mu = stats.t.pdf((t1 - my) / np.sqrt(s2/n), n-1) / np.sqrt(s2/n)\n",
    "\n",
    "# compute the exact marginal density for mu for the filtered data\n",
    "pm_mu_pos = (\n",
    "    stats.t.pdf((t1 - my_pos) / np.sqrt(s2_pos/n_pos), n_pos-1) /\n",
    "    np.sqrt(s2_pos/n_pos)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAI0CAYAAAAKi7MDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VPW9//HXhwQICRDBoCxBRIOgKLhQaq0KClIQN1yQ\npS5XLFytFm17UZQqYJVasFd/ar321q1aBRfcQGwtClivC1ZBUQpCixgCAoqACQlJ+P7+OCfDZMgy\nwdnOzPv5eORh5izf85lzHOaTz/f7Pcecc4iIiIgEWbNkByAiIiLyXSmhERERkcBTQiMiIiKBp4RG\nREREAk8JjYiIiASeEhoREREJPCU0ItIoM3Nm9uNkxyEiUh8lNCJNYGaP+l/uv41YXugvH5ik0JLO\nzNaZ2eXJjkOazszWmNnUZMch8l0ooRFpunLgZ2bWLdmBSOoys2ZmlpXsOEQyhRIakab7P2A5cEdD\nG5nZwX5FZ4uZ7TSzt8zs1LD1b5rZ7WGvp/lVnsFhy94ysxlhrwf7+5WZ2XYzW2xmh/vrzMx+aWb/\nMrPdZrbWzK6LiGmdmd1mZg+Y2TdmttnMrjGzlmZ2r5ltM7MNZnZNHW/pQDN7zsxK/W0mNvL+rzSz\nlWZWbmZfm9kSMytsYPtFZvaQmf3aj+sbM7vdTwxuMbMv/XN5e8R+zc1sqpn92z/WJ2Y2IWKbiWa2\nzMy+NbNNZjbbzDpFtPE7Mys2swoz22hms8PWP2pmf4to88dm5sJeT/UrHReb2T+B3cAR/rpR/vHL\n/WvwOzPLS9B7d2Z2tZk97v9/WGxmk8OPDRwO3Opv68zs0MbOiUjKcc7pRz/6ifIHeBT4G3AKsAfo\n5y8vBBww0H/dCvgUeA7oBxQBNwMVwJH+NtOBt8PafhPYDNzhv26N96V4hv96MFAN3A30BXoB44Be\n/vqfAruA8UAP4D/xqknjwo6xDvgG+Lkf0xQ/7lfClk3239tRYfs54GvgWrwv6YlAFXBuRNuX+7+f\n4K+/FOgGHANcCRQ2cG4XAduBO/1jXOEfdwHwW3/ZZf6yYRHX5CNgCNAduNh/j+Hve6J//roDP8BL\nSheHrf85UAwMBA4BvgdcF3ndI+L9sfdPaOj1VKAMWAx834+3DXA5sA24BDgMONWP9/EEvXcHfAn8\nBC9x+am/bJC/vj3wb2AW0NH/yWrsnOhHP6n2k/QA9KOfIP2Ef7EBzwOL/N8jE5rL/S+D7Ij9Xwfu\n9n8fCFT6X3q5eMnOL4B3/PXD/GWt/NdvAvMaiO0L4LcRy/4b+FfY63XAC2GvmwE7gJcjlm0Drglb\n5sK/gP1lTwJv1hPLCP8Lum0Tzu0iYFnEsk+AjyOWLQdm+b93x0u+ekVsc0tkWxHrj/PfUxf/9T3+\ntbHGrnvYsroSmj3AIRHbrQP+M2LZqf7x28X7vfvH+X8R26wEZoS9XgNMjdimwXOiH/2k2k82IrK/\nbgA+MbNzgA8i1n0P7y/db8wsfHlLvCoKwNt4VYwBeInN58DjwAwzawOcjpfc1Gx/AnBjXYGYWVu8\npGpJxKrFwEQzy3XOlfnLltesdM7tMbMteH/lhy/bDBwU0dbbEa/fAm6rKx7gNeBfwL/N7DW8L8a5\nzrmt9WxfY3nE603+T+Symtj6AQa8H3Ges/GqWQCYN1h7MnAUcAB7u9u7ARuAR/yY1/jxvoaX5O1u\nJN5IXzrn1ocdt4N/jN+Z2ayw7WqCLQKW+r/H5b37lkW8LgEObvCdxO6ciCSEEhqR/eScW21mD+J1\nEwyLWN0M76/gEXXsWubvX2Fm/wcMwutaet05t9nMVuElOacDL8Uh9MqI166eZfs9xs45962Z9QN+\niNfV85/Ab81skHPuHzGMrea/J+Gf14jtMLND8LrUHsfr5tuKl/z9DWjhx7vMzLoDZwCn4VUnbjOz\nE51zO/AqIRbRfvM64i+NeF0T30TgjTq2Lw77PebvPUxkEtLo9Y3inIikFA0KFvlupgGd8cathHsf\nb7zEDufcmoifkrDt3sBLXE4HFvrLXgcuAI71f6/xD7yxEvvwv2CK8boywg0A/h1WnfkuTox4fRLe\nOKE6OeeqnXNLnHO34FWXNgJjYhBHuJrk6JA6zvNaf9338MY0Xeece8s5t4o6qhPOuW+dc887536G\nV/04Eu/8gTe2qXPELsc3Fpxz7ku8rsCedcS3xjlX3uR3vFc07z1au/HGzUTG39A5EUkpqtCIfAfO\nuS1m9hvgVxGr/gxcD8w3s5uB1XhfoqcDK51zL/jbvY5XNahm71/wrwPP4g3ofSeszduABWZ2N/Aw\n3viaH+ANLF4FzADuMrPP8MZknA5chTcINBbO8mc//QUYijcA9aK6NjSzc/ESuiXAFryEpisNJED7\nwzm3xsweBv7XzCbhdYvl+cfr4Jy7E/gMryLxCzP7M96A6lsi4v0vvG6YZXjVjtF412S1v8nfgBvM\n7KfAq3jndmSUYd4MPGRm24AX8aouR+IN7p3Q4J4NiPK9R+vfwA/9alYZ3gDwX9DwORFJKarQiHx3\n/43XjRHi/+U9AK9S8wjel8BcoD/eWJkaS/G6KT4NG1+yGK974+/OucqwNv8KnIk3g+Zd4D28mS81\n2zyA90V9E17icANwo3PuoRi9z+l43UfL/WNMcs49X8+224Cz8b78V+PN1Pl1DGMJNx7vGtyM974X\n4p2XfwE45z7Cm501wV//S+C6iDZ24M3qeRv4GK+r8AI/UcQ59ze8GWE34b3/0/HOR6Occ4/jJT9n\n4V2zpXgDiDfsx3uN1OB7b4Jb8cYWrcJLQA+hkXMikmrMuciuVhEREZFgUYVGREREAk8JjYiIiASe\nEhoREREJPCU0IiIiEnhKaERERCTwAnsfmu3btwduelZubi5lZbG4v5nEm65VcOhaBYuuV3CkyrXK\nz8+PvEt3nVShSaCI561ICtO1Cg5dq2DR9QqOoF0rJTQiIiISeEpoREREJPCU0IiIiEjgKaERERGR\nwEtYQmNm15jZ+2ZWYWaPhi0/0cxeM7OvzWyLmT1jZp0SFZeIiIgEXyIrNCXAr4GHI5a3A/4AHAp0\nA3biPZ1YREREJCoJuw+Nc24ugJn1AwrDli8I387M7gMWJyouERERCb5UHENzKvBJsoMQERGR4Eip\nOwWbWR/gFuDcxrbNzc0N3E1/srKyyMvLS3YYEgVdq+DQtQqWIF2voil/iWq7Nb/+UZwjSY4gXStI\noYTGzIqABcBE59ybjW2fCrdjbqq8vDxKS0uTHYZEQdcqOHStgiUdr1e6vZ8aqXKt8vPzo9ouJbqc\nzKwb8DfgNufc48mOR0RERIIlYRUaM8v2j5cFZJlZDlAFHAy8DtznnPufRMUjIiIi6SORXU5TgFvD\nXv8YmAY44DBgqplNrVnpnGudwNhEREQkwBI5bXsqMLWe1dMSFYeIiIikn5QYQyMiIiLyXSihERER\nkcBTQiMiIiKBp4RGREREAk8JjYiIiASeEhoREREJPCU0IiIiEnhKaERERCTwlNCIiIhI4CmhERER\nkcBTQiMiIiKBl8iHU4qIiCRV3xl/T3YIEieq0IiIiEjgKaERERGRwFNCIyIiIoGnhEZEREQCTwmN\niIiIBJ4SGhEREQk8JTQiIiISeEpoREREJPCU0IiIiEjgKaERERGRwFNCIyIiIoGnhEZEREQCL2EJ\njZldY2bvm1mFmT0asW6Qmf3TzMrM7A0z65aouERERCT4ElmhKQF+DTwcvtDMCoC5wK+A9sD7wJwE\nxiUiIiIBl52oAznn5gKYWT+gMGzV+cAnzrln/PVTga1m1ss5989ExSciIiLBlQpjaHoDy2teOOdK\ngbX+chEREZFGJaxC04DWwJaIZduBNg3tlJubi5nFLah4yMrKIi8vL9lhSBR0rYJD1ypY4nG9iqb8\nJabtNVW6/v8XtM9WKiQ03wJtI5a1BXY2tFNZWVncAoqXvLw8SktLkx2GREHXKjh0rYIlHa9Xur2f\nGqlyrfLz86PaLhW6nD4B+ta8MLM84HB/uYiIiEijEjltO9vMcoAsIMvMcswsG3geONrMLvDX3wJ8\npAHBIiIiEq1EVmimALuAG4Ef+79Pcc5tAS4Abge2Ad8HRiUwLhEREQm4RE7bngpMrWfd34BeiYpF\nRERE0ksqjKERERER+U6U0IiIiEjgKaERERGRwFNCIyIiIoGnhEZEREQCTwmNiIiIBJ4SGhEREQk8\nJTQiIiISeEpoREREJPCU0IiIiEjgKaERERGRwFNCIyIiIoGnhEZEREQCTwmNiIiIBJ4SGhEREQk8\nJTQiIiISeEpoREREJPCU0IiIiEjgKaERERGRwFNCIyIiIoGnhEZEREQCTwmNiIiIBJ4SGhEREQk8\nJTQiIiISeCmT0JjZoWb2ipltM7NNZnafmWUnOy4RERFJfSmT0AC/BzYDnYBjgQHA1UmNSERERAIh\nlRKa7sDTzrly59wm4FWgd5JjEhERkQBIpYTmbmCUmeWaWRdgGF5SIyIiItKgVBqjsgQYD+wAsoDH\ngBfq2zg3NxczS1BosZGVlUVeXl6yw5Ao6FoFh65VsKTj9Uq391MjaNcqJRIaM2uGV435A3AS0Bp4\nGLgTmFTXPmVlZQmLL1by8vIoLS1NdhgSBV2r4NC1CpZ0vF7p9n5qpMq1ys/Pj2q7VOlyag8cAtzn\nnKtwzn0FPAKcmdywREREJAhSIqFxzm0F/g1cZWbZZnYAcBnwUXIjExERkSBIiYTGdz4wFNgCrAEq\ngeuTGpGIiIgEQkqMoQFwzi0DBiY7DhEREQmeVKrQiIiIiOwXJTQiIiISeEpoREREJPCU0IiIiEjg\nKaERERGRwFNCIyIiIoGnhEZEREQCTwmNiIiIBF7UCY2Z/bKe5T+PXTgiIiIiTdeUCs0t9SyfEotA\nRERERPZXo48+MLPT/V+zzOw0wMJWHwbsjEdgIiIiItGK5llOD/n/zQEeDlvugE3AtbEOSkREJB31\nnfH3qLddPvnkOEaSfhpNaJxz3QHM7E/OuUvjH5KIiIhI00T9tO3wZMbMmkWs2xPLoERERESaoimz\nnI43s7fNrBSo9H+q/P+KiIiIJE3UFRrgMeBl4AqgLD7hiIiIiDRdUxKabsDNzjkXr2BERERE9kdT\n7kPzPDAkXoGIiIiI7K+mVGhygOfN7O9407VDNPtJREREkqkpCc2n/o+IiIhISmnKtO1p8QxERERE\nZH9FndCEPQJhH86512MTjoiIiEjTNaXL6aGI1x2AFkAx3jOdRERERJKiKV1O3cNfm1kW3pO29XBK\nERERSaqmTNuuxTlXDdwOTIpVMGY2ysxWmlmpma01s1Ni1baIiIikr6Z0OdXlDCAmz3EyszOAO4GL\ngfeATrFoV0RERNJfUwYFfwGE3yU4F+/eNFfHKJZpwHTn3Dv+6w0xaldERETSXFMqND+OeF0KrHbO\n7fiuQfjjcfoBL5nZGrxE6QXgv5xzu75r+yIiIpLemjIoeDGAmTUDDga+dM7FpLvJb685cCFwCt4T\nvF/EG3R8c1075ObmYmYxOnxiZGVlkZeXl+wwJAq6VsGhaxUs6Xi94vV+kn2egnatmtLl1Aa4H2+M\nS3Og0sxmAz9zzm3/jnHUVGHudc5t9I/3OxpIaMrKgvfA77y8PEpLS5MdhkRB1yo4dK2CJR2vV7ze\nT7LPU6pcq/z8/Ki2a8osp3uBPOAYoJX/31zg/zU1uEjOuW1497MJH6Ojp3qLiIhIVJoyhmYocJhz\nrqY0strM/gNYG6NYHgGuNbNX8bqcrgfmxahtERERSWNNqdCU490dOFwBUBGjWG4DlgKrgZXAh3j3\nuRERERFpUFMqNH8EXvPHtnwOdMOrovxvLAJxzlXiTQGP1TRwERERyRBNSWhux7s3zFigM1AC/NY5\nF/mMJxEREZGEakqX0z3AKufcYOfcUc65wcBKM7s7TrGJiIiIRKUpCc1o4P2IZf8AxsQuHBEREZGm\na0pC44CsiGVZTWxDREREJOaakoy8Cdzm3ym45o7BU/3lIiIiIknTlEHBE/HuC7PRzD4HDgE2AmfH\nIzARERGRaDXlWU7FZnY80B/oCnwBvBfD5zmJiIiI7JemVGjwk5d3/B8RERGRlKABvSIiIhJ4SmhE\nREQk8JTQiIiISOA1aQyNiIhIIvSd8fdkhxC1VIg12hiWTz45zpEkjyo0IiIiEnhKaERERCTwlNCI\niIhI4CmhERERkcBTQiMiIiKBp4RGREREAk8JjYiIiASeEhoREREJPCU0IiIiEnhKaERERCTwlNCI\niIhI4CmhERERkcBLuYTGzHqYWbmZPZHsWERERCQYUi6hAe4HliY7CBEREQmOlEpozGwU8A2wMNmx\niIiISHBkJzuAGmbWFpgOnA5c2dj2ubm5mFnc44qlrKws8vLykh2GREHXKjh0rYJF1yt68ThPTWkz\naNcqZRIa4DbgIedccTSJSllZWfwjirG8vDxKS0uTHYZEQdcqOHStgkXXK3rxOE9NaTNVrlV+fn5U\n26VEQmNmxwKDgeOSHYuIiIgET0okNMBA4FBgvV+daQ1kmdlRzrnjkxiXiIiIBECqJDR/AGaHvf4l\nXoJzVVKiERERkUBJiYTGOVcGhAbFmNm3QLlzbkvyohIREZGgSImEJpJzbmqyYxAREZHgSKn70IiI\niIjsDyU0IiIiEnhKaERERCTwlNCIiIhI4CmhERERkcBTQiMiIiKBp4RGREREAk8JjYiIiASeEhoR\nEREJPCU0IiIiEngp+egDERFJT0VT/pLsECRNqUIjIiIigaeERkRERAJPCY2IiIgEnhIaERERCTwl\nNCIiIhJ4SmhEREQk8JTQiIiISOApoREREZHAU0IjIiIigaeERkRERAJPCY2IiIgEnhIaERERCTwl\nNCIiIhJ4KZHQmFlLM3vIzD43s51mtszMhiU7LhEREQmGlEhogGzgC2AAkA9MAZ42s0OTGJOIiIgE\nRHayAwBwzpUCU8MWzTOzfwMnAOuSEZOIiIgER0okNJHM7GDgCOCT+rbJzc3FzBIXVAxkZWWRl5eX\n7DAkCrpWwaFrFR9FU/4S9bZrfv2jOEaSueLx/3VT2gzaZyvlEhozaw78GXjMOffP+rYrKytLXFAx\nkpeXR2lpabLDkCjoWgWHrlXy6fzHRzzOa1PaTJXPVn5+flTbpcoYGgDMrBnwOLAbuCbJ4YiIiEhA\npEyFxrz+o4eAg4EznXOVSQ5JREREAiJlEhrgAeBIYLBzbleygxEREZHgSIkuJzPrBkwAjgU2mdm3\n/s/YJIcmIiIiAZASFRrn3OdAsKYsiYiISMpIiQqNiIiIyHehhEZEREQCTwmNiIiIBJ4SGhEREQk8\nJTQiIiISeEpoREREJPCU0IiIiEjgKaERERGRwFNCIyIiIoGnhEZEREQCz5xzyY5hv2zfvj2ugfed\n8feotls++eSo28zLy6O0tHR/Q5IE0rVKvmg/g2t+/SNdqziI9vyLQNO+C5sqPz8/qkcjqUIjIiIi\ngaeERkRERAJPCY2IiIgEnhIaERERCTwlNCIiIhJ4SmhEREQk8JTQiIiISOApoREREZHAU0IjIiIi\ngaeERkRERAJPCY2IiIgEnhIaERERCTwlNCIiIhJ4KZPQmFl7M3vezErN7HMzG5PsmERERCQYspMd\nQJj7gd3AwcCxwHwzW+6c+yS5YYmIiEiqS4kKjZnlARcAv3LOfeuc+zvwEnBJciMTERGRIDDnXLJj\nwMyOA95yzuWGLfslMMA5d3byIhMREZEgSIkKDdAa2BGxbDvQJgmxiIiISMCkSkLzLdA2YllbYGcS\nYhEREZGASZWEZjWQbWY9wpb1BTQgWERERBqVEmNoAMxsNuCAK/FmOb0CnKRZTiIiItKYVKnQAFwN\ntAI2A08BVymZERERkWikTIVGREREZH+lUoVGREREZL8ooREREZHAU0KTQGbWw8zKzeyJiOVj/OdX\nlZrZC2bWPlkxZjIza2lmD/nXYqeZLTOzYRHbDDKzf5pZmZm9YWbdkhVvptPz31JXY58lfY5SU13f\nUUH6flJCk1j3A0vDF5hZb+BBvMc8HAyUAb9PfGiC92yzL4ABQD4wBXjazA4FMLMCYC7wK6A98D4w\nJxmBClD7+W9jgQf8z5MkX72fJX2OUlqt76igfT9pUHCCmNko4HzgU6DIOfdjf/kdwKHOuTH+68OB\nlcCBzjndWDDJzOwjYJpz7jkzGw9c7pw7yV+XB2wFjnPO/TOZcWYa/9xvA452zq32lz0ObHDO3ZjU\n4KRONZ8l4ED0OUo5dX1HBe37SRWaBDCztsB04Od1rO4NLK954Zxbi/dX5xGJiU7qY2YH412HmtsH\nRF6rUmCtv1wS6wigqiaZ8S1H1yIlRXyW9DlKMQ18RwXq+0kJTWLcBjzknCuuY11rvOdWhdNzrJLM\nzJoDfwYeC/urUdcqdej5bwFRx2dJn6PUU993VKCulRKa78jMFpmZq+fn72Z2LDAY+O96mtBzrBKk\nsWsVtl0z4HG8v0SuCWtC1yp16FoEQD2fJV27FNLId1SgrlV2sgMIOufcwIbWm9l1wKHAejMDL+PN\nMrOjnHPH45Vg+4ZtfxjQEu/5VhJDjV0rAPMu0kN4A+DOdM5Vhq3+BLgsbNs84HD0zLFkCD3/zTn3\nmb9Mz39LIQ18lvQ5Si0Dqec7CniVAH0/aVBwnJlZLrUz3F/i/c9zlXNuiz+K/G1gOPAB3ojybOfc\nqETHKmBm/4P3LLHBzrlvI9Z1ANYAVwDz8QY4DnDOnZjwQEXPf0tx9X2W9DlKLQ19RwEHEaDvJ3U5\nxZlzrsw5t6nmB6+EV+6c2+Kv/wT4T7w+5s14fZNXJy3gDObfC2MC3j/Cm8zsW/9nLIB/zS4Abseb\nYfN9ICU/2BlCz39LUQ19lvQ5Si0NfUcF7ftJFRoREREJPFVoREREJPCU0IiIiEjgKaERERGRwFNC\nIyIiIoGnhEZEREQCTwmNiIiIBJ4SGhEREQk8JTQiIiISeEpoRCSlmdnN/m30a163M7NKM8tJZlwi\nklqU0IhIqjsGWBb2+lhglXOuPEnxiEgKUkIjIqmuroRmeZJiEZEUpYRGRFKWmbUADgc+Clvcl9oJ\njoiIEhoRSWlHAhucc2UAZmbAQFShEZEISmhEJJX1AQ4ys8PNrBVwG9ANWJfUqEQk5SihEZFUdgzw\nF2ARsAbYCRQDNycxJhFJQeacS3YMIiJ1MrMFwB+dc88lOxYRSW2q0IhIKjsGWJnsIEQk9alCIyIp\nyczaAV8Cec65ymTHIyKpTQmNiIiIBJ66nERERCTwlNCIiIhI4CmhERERkcBTQiMiIiKBp4RGRERE\nAk8JjYiIiASeEhoREREJPCU0IiIiEnhKaERERCTwlNCIiIhI4CmhERERkcBTQiMiIiKBp4RGRERE\nAk8JjYiIiASeEhoRSXlmtsjM/pjsOGqY2UAzc2ZW2IR9LjezqnjGJZLJlNCIZAAze9T/Av5txPJC\nf/nAJIUmIhITSmhEMkc58DMz6xbLRs3TPJZtiog0lRIakczxf8By4I6GNjKznmY238y+9X9eNrOi\nsPWXm1mVmZ1mZh8CFcBgM5tqZmvMbKSZfWZmZWb2gpm1NbPzzWyVme00s2fNLD+svePNbIGZbfaP\nt9TMhjbljZnZoX6laYyZ/cU/9j/NbICZdTGzV8ys1Mw+NbNTIvY90cyWmNkuM9tmZk+a2UER21xr\nZsV+u38BDqkjhhPM7K/+e9hiZnNjnTyKSP2U0IhkDgf8EhhtZv3q2sDMWgF/BXKAAf5Pa+BVM2sR\ntmkz4E7g50Av4H1/eSfgMuACYBjwQ+BZ4EpgpL/sFOCmsLbaAnOA04Djgb8AL5nZEfvxHm8DHgCO\nBVYCs4HHgP8FjgM+BZ6sqSiZWUf//RYD/YGzgaP9mGvOybnAfwO/89t9GpgZflAzOwpYDLwN9ANO\nB6qB18wsZz/eh4g0UXayAxCRxHHOvWlmLwKzgIF1bDIG6ACc4JzbCmBmo4B1wCjgT/52BvzCOfdm\nzY5mBtASuCxs36eB/wQ6Oue2+MtmA4PCYloUEcMUMzsbuAi4vYlv8V7n3Av+ce4A3gNmOeee95fd\nDnwA9ARWAD8FdgCXO+d2+9tcAiwzs1Odc0uA/wLmOOd+5x9jtZkdCfwi7LiTgHnOuVvDzsePgW3A\nUOCFJr4PEWkiVWhEMs8NwA/N7Jw61vUGPq1JSACcc18Cq/x14ZbWsf+G8H2BTcCmmmQmbFmoS8fM\nOpjZ7/0uom/M7Fv/WPvTXbM84jgAH9WxrOb4vYF3apIZAOfccmA7e9/vUXjddeH+HvH6e8CIsG66\nb4Gv8CpdPfbjfYhIE6lCI5JhnHOrzexBvC6jYfvZTLVzrryO5ZWRh6tnWfgfU4/ijUmZBPwb2IXX\nVdSCpgs/lmtgWaz/mGsGPA78po51X8X4WCJSB1VoRDLTNKAzMD5i+SfAUWZWULPAzA5mbxdNPJwK\n/N4595Jz7mNgI3BYnI4V6RPgxPDxQWbWF8hn7/v9FDgpYr8fRrx+H+gDrHXOrYn42Ran2EUkjBIa\nkQzkdwH9BrguYtWTwBZgjj/76AS8askGvIG78bAKGGtmx5jZscBTQFacjhXpPrxByY+a2dFmdjJe\npeXNsPFBdwEXm9lEM+thZv8BXBLRzh3AkcATZtbfzLr7s8DuMbNEJWciGU0JjUjm+m8gfLwLzrld\nwBC8qdhL8GbulAJDw8eZxNh/4P1b9B7e4NlXqXt8Tsz544OGAIX+MefhVWYuDNvmebwBwJPwxuOM\nxRuHFN7OSrwqTmu8WVqf4s2sagV8E+/3ISJgzrnGtxIRERFJYarQiIiISOApoREREZHAU0IjIiIi\ngaeERkRERAIvsDfW2759e+BGM+fm5lJWVpbsMCQKulbBoWsVLLpewZEq1yo/P9+i2U4VmgTyn3Uj\nAaBrFRy6VsGi61XbRx99xEcffdT4hkkQtGsV2AqNiIhI0E2ePBmA+fPnJzmS4FOFRkRERAJPCY2I\niIgEnhIaERERCby0HENTVlbGjh07qKqqSnYotZgZetREMAT9WmVnZ9O2bVtyc3OTHYqISEKkXUJT\nVlbG9u2xodp/AAAgAElEQVTb6dKlCzk5OSk1SjsrK4vq6upkhyFRCPK1cs5RXl7Ohg0bAJTUiKSw\nW265JdkhpI20S2h27NhBly5daNWqVbJDEUkKM6NVq1Z06dKFkpISJTQiKez73/9+skNIG2k3hqaq\nqoqcnJxkhyGSdDk5OSnX7SoSBLajGNtZkpBjvfvuu7z77rsJOVa6S7sKDQTvZkAi8aDPgch+2FNF\n62cuprp9EWUjHov74aZPnw7oPjSxkHYVGhERkf3VfPU8mu34AivfluxQpImU0IiIiAA4R8v3/wcA\nq9yV5GCkqZTQiIiIAFlf/B9ZWz7FtWiNVSb/oYzSNEpoJKRv374sXrw47sdZtWoVJ5xwAu3atePe\ne++N+/FERKKRtW0tAFWFP4AqVWiCJi0HBWeioqIiHnzwQQYNGrTfbSxfvjyGEdVv1qxZDBw4kH/8\n4x8JOZ6ISFSqygFwuQcmrMtpxowZCTlOJlCFRr7T1N792Xf9+vUcddRR+31MEZF4qEli9rRqD5Vl\nkIC7hffp04c+ffrE/TiZQAlNEhQVFXHnnXfSp08fOnTowLhx4ygv9/4yWLlyJYMGDaKgoIC+ffvy\n8ssvh/abOXMm3bp1o127dvTu3ZvXX38dgMsuu4z169dz3nnnccABBzBr1ixKSkoYOXIknTp1okeP\nHvt07RQVFTFz5kyOO+448vPzqaqqoqioiIULFzYaR137Rqpv/zPOOINFixYxceJEDjjgAFavXl3n\nOZoxYwZXX3116PW2bdto1apV6DyJiMRcVTmuWTa0aIPhoLoi7odctGgRixYtivtxMkHadzn9ev5K\nPt24M67HOKpTG6YMP7JJ+zz11FPMnz+fvLw8zjvvPO644w5+9atfMWLECC6//HIWLFjAW2+9xfnn\nn88777wDwO9//3vefvttOnfuzLp160K35n/sscd46623Ql1Oe/bs4cQTT+Scc87hiSeeoLi4mKFD\nh9KzZ0+GDBkSimH27Nm89NJLFBQUkJ2993+FysrKeuPo2bNng/s2tv9rr73GoEGDGDNmDOPGjav3\n/KxYsYJTTz019Hr58uX07NlTN00Ukbixql2QnYNr7t1p3ip34bLj+2/OzJkzARg4cGBcj5MJVKFJ\nkquuuoquXbvSvn17Jk+ezJw5c3j33Xf59ttvmTRpEi1atOC0005j+PDhzJkzh6ysLCoqKli5ciWV\nlZUceuihHH744XW2vXTpUrZu3cqUKVNo0aIFhx12GOPGjWPOnDm1trvmmmvo2rXrPo+JaCiOxvaN\ndv/GrFixgr59+4ZeL1u2jGOOOSbq/UVEmqyqHJfdCtfcf1yIZjoFStpXaJpaOUmUrl27hn4/5JBD\nKCkpoaSkhMLCQpo1a1Zr3YYNGygqKuKuu+5i+vTpfPrppwwZMoSZM2fSuXPnfdpev349JSUlFBQU\nhJZVV1dz8skn1xtDuIbiaGzfaPdvyO7du1m7dm2tfuWPPvqIY489Nqr9RUT2h1V6FRr8hMYqy4j/\nKBqJFVVokuSLL76o9Xvnzp3p3LkzxcXF7Nmzp9a6Ll26ADB69GgWL17M2rVrMTMmT54c2i78NveF\nhYV0796drVu3hn62bdtWaxxM5D7hGoujoX2j3b8hK1eupEuXLqGHKjrnWLJkiQbOiUhcWVU5Ljtn\nbzeTpm4HihKaJHnggQcoLi7m66+/ZsaMGVx00UX079+f3NxcZs2aRWVlJYsXL2bevHmMHDmSVatW\n8cYbb1BRUUFOTg45OTm1KiAHHXQQ//rXvwDo378/rVu3ZubMmezatYvq6mpWrFjB0qVLo4qtoTgS\nsf/HH3/M5s2bWbt2Lbt27eLWW2/l888/p1u3blHtLyKyX2oSmlCFRglNkCihSZLRo0dz5plncsQR\nR3DYYYdx00030aJFC55//nleffVVOnbsyLXXXssjjzxCr169qKio4KabbqJjx44UFhayZcsWbr/9\n9lB7N9xwAzNmzKCgoIB77rmHF198keXLl9OjRw86duzIhAkT2LFjR1SxNRRHIvZfsWIFQ4YMYfDg\nwfTq1Ys2bdpQWFio+zWISFxZVTk0b1Wryyne7r77bu6+++64HycTmEvAPPt42L59e52BFxcXh2bi\npJqsrCyqq6tjchO8dHbWWWdxxRVXcP755ycthpprFXSrVq2isLAw2WHEVV5eHqWlpckOQ6KUytcr\n76lzcC3bUn7KzbR54keUnvU/VPU4M9lhJU2qXKv8/Pz6xziEUYVGUs6KFSuiruaIiMSKVZVDdiuv\nSkNiKjQLFixgwYIFcT9OJkj7WU4SLNu2bWPz5s306NEj2aGISKbx7zuTyDE09913HwDDhg2L+7HS\nnRKaJFizZk2yQ0hZ7dq1o6xM934QkcSrGUPjsv37a2mWU6Coy0lERATvTsEu4j40EhxKaERERACq\nKrwb6zXLwmW1VEITMEpoRERE9lRj1RWh7iaXnQNVehhukGgMjYiIiP9k7dBdgpvnJqRC8+CDD8b9\nGJlCCY2IiGS80Iwmf8q2a56bkIdTpvt9ohJJXU4iIiJ+99LeCk2rhEzbnjt3LnPnzo37cTKBKjQi\nIpLxrGaKdmgMTau9y+LooYceAkjqndHThSo0IiIiERWaRHU5SewooRERkYwXOYYmUV1OEjtKaCSk\nb9++LF68OO7HWbVqFSeccALt2rXj3nvvjfvxEvW+RCTAaio0WS29/2a30n1oAkYJTZooKipi4cKF\n36mN5cuXM2DAgBhFVL9Zs2YxcOBAtm3bxrXXXhv34yXqfSXC119/zYUXXkh+fj6HH344Tz31VLJD\nEkkLNeNlXPgsJz36IFA0KFioqqoiO3v//lfYn33Xr1/PyJEj9+t46eTLL7/k4IMPbtI+P/vZz2jR\nogUbNmxg2bJlnHvuufTp04fevXvHKUqRzGA1N9FL8CynP/3pT3E/RqZQhSYJioqKuPPOO+nTpw8d\nOnRg3LhxlJd7H6aVK1cyaNAgCgoK6Nu3Ly+//HJov5kzZ9KtWzfatWtH7969ef311wG47LLLWL9+\nPeeddx4HHHAAs2bNoqSkhJEjR9KpUyd69OixT9dOUVERM2fO5LjjjiM/P5+qqqpaVZ6G4qhr30j1\n7X/GGWewaNEiJk6cyAEHHMDq1avrPEczZszg6quvDr3etm0brVq1Cp2ncB988AH9+vWjXbt2jBo1\nijFjxnDLLbfUinfhwoXMnDmTiy++uNa+119/Pddddx1ArXN22GGH1TpnRUVF/O53v+O4447jwAMP\nZMyYMXXGAvDtt9/SsmVLNm7cGFq2YsUKunbtys6dO0PLxo0bxw9+8AMefPBBvvnmmzrbCldaWsrc\nuXOZOnUqrVu35uSTT+ass87iz3/+c6P7ikgjQoOCI+5D41xcD3vggQdy4IEHxvUYmSLtKzTN/noT\ntmlFXI/hOh7NniF3NGmfp556ivnz55OXl8d5553HHXfcwa9+9StGjBjB5ZdfzoIFC3jrrbc4//zz\neeeddwD4/e9/z9tvv03nzp1Zt24d1dXVADz22GO89dZbPPjggwwaNIg9e/Zw4okncs455/DEE09Q\nXFzM0KFD6dmzJ0OGDAnFMHv2bF566SUKCgpqVVkqKyvrjaNnz54N7tvY/q+99hqDBg1izJgxjBs3\nrt7zs2LFCk499dTQ6+XLl9OzZ09ycnJqbbd7924uuugiJk6cyFVXXcW8efMYO3Ysv/zlL/dpc+TI\nkdx2223s3LmTNm3aUF1dzbPPPsszzzzDnj17OO+880LnbOPGjZxxxhm1ztkzzzzD/PnzycnJYcCA\nATz22GNMmDBhn+O0bt2aXr168eGHH9KpUycAbr75Zm644QbatGkT2u7555/nlVde4fHHH+emm27i\nzDPP5LLLLuP000+nWbN9/9ZYvXo12dnZHHHEEaFlffv2ZcmSJfWeRxGJzj4VmuxWGM67g3B2Tv07\nfkc1f5CMHTs2bsfIFKrQJMlVV11F165dad++PZMnT2bOnDm8++67fPvtt0yaNIkWLVpw2mmnMXz4\ncObMmUNWVhYVFRWsXLmSyspKDj30UA4//PA62166dClbt25lypQptGjRgsMOO4xx48YxZ86cWttd\nc801dO3alVatWtVa3lAcje0b7f6NWbFiBX379g29XrZsGcccc0ydx6qqquLaa6+lefPmjBgxgu99\n73t1ttmtWzeOO+44XnjhBQDeeOMNcnNzOfHEE6M6Z9dccw2dO3emffv2DB8+nOXLl9cb/wknnMCH\nH34IwJtvvsnKlSsZP358rW2aN2/Oueeey7PPPsuqVavo378/kydPpqioiPvvv3+fNktLS2nbtm2t\nZW3btq1V9RGR/VQZOYbG+2+8u52efPJJnnzyybgeI1OkfYWmqZWTROnatWvo90MOOYSSkhJKSkoo\nLCys9df5IYccwoYNGygqKuKuu+5i+vTpfPrppwwZMoSZM2fSuXPnfdpev349JSUlFBQUhJZVV1dz\n8skn1xtDuIbiaGzfaPdvyO7du1m7di19+vQJLfvoo4849thj6zxW586dMbOoYhs1ahRz5szhkksu\n4amnnmLUqFFAdOesY8eOod9zc3MpKSmp9zj9+vULdQneeOONTJ06lRYtWtS7/YEHHkifPn3o27cv\nzz33HOvWrdtnm7y8PHbs2FFrWU21SUS+m8gKjWue672uLINW7ZIUlTSFKjRJ8sUXX9T6vXPnznTu\n3Jni4mL27NlTa12XLl0AGD16NIsXL2bt2rWYGZMnTw5tF/6FXlhYSPfu3dm6dWvoZ9u2bbXGwUTu\nE66xOBraN9r9G7Jy5Uq6dOlCbq73D4pzjiVLltRKcGp06tSJkpISXFg/d/i5jXThhReyePFiiouL\nefHFF0MJTeQ527ZtW53nLFr9+vXjww8/ZO7cuZSXlzN69Og6t/vss8+49dZb6dGjB9dffz1HH300\nq1evZubMmftse8QRR1BVVcVnn30WWrZ8+XKOOuqo/YpRRMJU7cI1aw7N/L/z/YRGU7eDQwlNkjzw\nwAMUFxfz9ddfM2PGDC666CL69+9Pbm4us2bNorKyksWLFzNv3jxGjhzJqlWreOONN6ioqCAnJ4ec\nnJxaFZCDDjqIf/3rXwD079+f1q1bM3PmTHbt2kV1dTUrVqxg6dKlUcXWUByJ2P/jjz9m8+bNrF27\nll27dnHrrbfy+eef061bt322PfHEE8nKyuL++++nqqqKl156qcH32aFDBwYMGMCVV17JoYceypFH\nHhmK+bucs0h9+/Zl06ZNTJo0idtvv73OBPDKK6/klFNO4ZtvvuHpp5/mgw8+4LrrrqNDhw51tpmX\nl8eIESOYNm0apaWlvPXWW7z88svqexeJAasqrzVWJvRMJ03dDgwlNEkyevRozjzzTI444ggOO+ww\nbrrpJlq0aMHzzz/Pq6++SseOHbn22mt55JFH6NWrFxUVFdx000107NiRwsJCtmzZwu233x5q74Yb\nbmDGjBkUFBRwzz338OKLL7J8+XJ69OhBx44dmTBhwj7dFfVpKI5E7L9ixQqGDBnC4MGD6dWrF23a\ntKGwsJAZM2bUeaynn36aRx55hIKCAp588kmGDx9Oy5Yt621/1KhRLFy4sFbVJCsrq9Y569ChQ5PO\nWaSWLVty9NFH061bN4YOHVrnNuPHj2f9+vXcc889HH/88VG1e++997Jr1y46d+7MJZdcwn333acp\n2yIxYJW79iYx7O1y0t2Cg8NcnKekxcv27dvrDLy4uDg0EyfVZGVlUV1dTVFRUWhGkuzrrLPO4oor\nrtjvh7WddNJJjB8/nssvv3y/Y6i5Vvtr9+7d9OrViyeffJITTzxxv9v5rlatWkVhYWHSjp8IeXl5\nlJaWJjsMiVKqXq9WCyaSXfI+O8e9BUDWxg9oPfs8Ss97jKrup8XtuGVlXpdWTRd7KkmVa5Wfn1//\nGIcwqtBIylmxYkXU1RyAJUuWsGnTJqqqqvjTn/7Exx9/zI9+9KM4Rti42267jZNOOimpyYyIRM+q\nykMzm2Dv/Wji3eWUm5ubkslMEKX9LCcJlm3btrF582Z69OgR9T6rVq1i9OjRlJaW0r17d+bMmRO6\n/0uiffDBB5xxxhkcc8wxPPvss0mJQUT2Q1V5rS4nQtO24zso+I9//CPgjamT70YJTRKsWbMm2SGk\nrHbt2oVKsNH6yU9+wk9+8pM4RdQ0xx9/PF999VWywxCRJrKqXZAdVqHxH1JJ9e64Hvf5558HlNDE\ngrqcREREIis0Wd59oyzOCY3EjhIaERHJeFa5K9TNBOD8hCbeFRqJHSU0IiIiqtAEnhIaERHJeN4Y\nmn0TGqoqkhOQNJkGBYuISMazyAqNGS6rBVYd34Rm/vz5cW0/k6hCIyIiUlW+994zNbJaaAxNgCih\nERGRzLan2hsrE16hwZ+6HeeE5t577+Xee++N6zEyhRIaaZJ169bRvHlzqqqqkh2KiEhs+EmLy454\nBlxWi7gPCn711Vd59dVX43qMTKGEJsGKiopYt24dV1xxBY899hgAjz32GAMGDEhyZIk3aNAgFi9e\nzPTp05k+fXqywxGRTFUzTqZmILDPqcspUJTQBMR3eVCiiIjUL1SFyaqjQqNZToGhhCbJVq5cyU9/\n+lPeeecdDjjgAAoKCgC44oor+OlPf8rZZ59Nfn4+ixYtYtCgQTz00EOhfSMrO//85z8ZOnQoBx10\nEL179+aZZ56p85hPP/003//+92stu/vuuxkxYgQAr7zyCv369aN9+/Z07969wepJUVERCxcuDL2e\nPn06l156aej1O++8wymnnEJBQQHHH388ixcvbsLZERFJgJoup4gKDdnxH0MjsZMR07YHDRq0z7IL\nL7yQq666irKyMs4+++x91l966aVcdtllbN26lYsvvnif9RMmTGDkyJF88cUXdO3aNepYap7j9PDD\nD4eW3X///Tz88MP7fNnPnj2bl156iRdffJHduxv+UJWWljJs2DBuvfVW5s2bx8cff8ywYcPo3bs3\nRx11VK1tzzrrLCZMmMBnn30Wegjk7Nmzuf766wHvkfGPPPIIvXv3ZsWKFQwbNoy+ffty7rnnRv0+\nATZs2MC5557Lo48+yo9+9CNef/11Ro4cyYoVK+jQoUMoEcrE7jYRSR17KzSJ73Jq1apV4xtJVFSh\nSWFnn302P/zhD2nWrBk5OTkNbjt//ny6devG5ZdfTnZ2NscddxwjRozgueee22fb3Nxczj77bObM\nmQPAZ599xqpVq0KJ3YABAzjmmGNo1qwZffr04eKLL2bJkiVNjv/JJ59k6NChDBs2jGbNmjF48GBO\nOOEEFixY0OS2RETipr4KTQLuQ/Pss8/y7LPPxvUYmSIjKjThXSKRcnNzG1xfUFDQ4PqmVGeaqilt\nr1+/nvfeey/UZQVQVVXF2LFj69x+9OjRTJo0iSlTpjB79mzOOecccnNzAXj33Xe5+eab+eSTT9i9\nezcVFRVccMEFTY7/888/57nnnqt146jKykoGDhzY5LZEROIlNE6mjgqNVexIQkSyPzIioUl1ZhbV\n8ry8PMrKykKvv/zyy9DvhYWFnHrqqVFP/xs8eDBbtmxh2bJlzJkzh1mzZoXWXXrppVx11VXMmzeP\nnJwcfv7zn7N169Y624mMadOmTaHfu3btytixY3nwwQejiklEJClqqjD7TNtuGfdp27/97W8BmDRp\nUlyPkwnU5ZQCDj74YDZs2NDoOJm+ffvywgsvUFZWxpo1a3jkkUdC64YPH85nn33GE088QWVlJZWV\nlSxdupSVK1fW2Vbz5s254IILuPHGG/n6668ZPHhwaN3OnTtp3749OTk5vPfee8yePbvemPr06cPT\nTz9NZWUl77//PnPnzg2tGzNmDPPnz+evf/0r1dXVlJeXs3jxYoqLi6M9NSIi8VdPl1MixtAsXrxY\nkyViRAlNCjjttNM46qijKCwspGPHjvVuN3HiRFq0aEGXLl244oorGD16dGhdmzZteOWVV3j66ac5\n5JBDKCws5KabbqKiov7+39GjR7Nw4UIuuOACsrP3Fuvuvfdepk2bRrt27bj99tu58MIL621j2rRp\nrF27lg4dOjB9+nRGjRoVWte1a1eee+45fvOb39CpUye6d+/OXXfdxZ49e6I9NSIicadp2+nBnHPJ\njmG/bN++vc7Ai4uL6dmzZ6LDiUpWVpbuJxMQ6XKtVq1aRWFhYbLDiKu8vDxKS0uTHYZEKRWvV/aa\nV8l7eTw7xy5gz0G9Q8tbvTaJ7H+/wc7xS+N27OHDhwOp+ZDKVLlW+fn5dY/LiKAKjYiIZLRkTtuW\n2NGgYBERyWx+t1Ld07bjm9C0b98+ru1nEiU0IiKS2WqSlohZTl6FJr5jaB5//PG4tp9J0rLLKajj\ngkRiSZ8DkejU1+VEVktsTxU4TWQIgrRLaLKzsykvL092GCJJV15eXmv2mojUo4E7BYevj4dp06Yx\nbdq0uLWfSdLuX7u2bduyYcMGunTpQk5OTr03rRNJV845ysvL2bBhA/n5+ckORyTl1TsouKYLqqoC\nsht+/Mz+eu+99+LSbiZKu4Sm5vb9JSUlVFVVJTma2sxM3QABEfRrlZ2dTX5+fujzICINqElomjWv\nvdxPcKx6N8H91yBzpF1CA15Sk4r/kKfKnH5pnK6VSOawqnJcVkuIqOi7BHQ5Seyk3RgaERGRJqne\nve+AYKhVoZHUl5YVGhERkahV7953QDDhFZr4Td3u0qVL3NrONEpoREQko1n17n2ftA17KzRxfJ7T\nH/7wh7i1nWnU5SQiIpmtngpNaGaTupwCQQmNiIhkNKtnDE0iBgXfeOON3HjjjXFrP5Ooy0lERDJb\ndUXSBgV//PHHcWs706hCIyIima16tzdtO4KmbQeLEhoREclo9XU57a3QxPcBlRIbSmhERCSz1Tco\nWBWaQNEYGhERyWhW731owp7lFCdFRUVxazvTKKEREZHM1kiFJp6Dgu+55564tZ1p1OUkIiIZzarq\nnuUUetq2upwCQQmNiIhktuqKhu8UHMeEZuLEiUycODFu7WcSdTmJiEhmq2faNs2ah9bHy5o1a+LW\ndqZRhUZERDJavdO2zXBZLTVtOyCU0IiISGarb1AweImOxtAEghIaERHJXG4Ptqeq7goN/t2CldAE\ngsbQiIhI5qpJVhqo0Fgc70NzzDHHxK3tTKOERkREMpefrLi6ZjkR/wrNb37zm7i1nWnU5SQiIhnL\nGqvQZLeM67RtiR0lNCIikrn8ZKXhQcHx63IaP34848ePj1v7mURdTiIikrEaq9C4rJZx7XLasGFD\n3NrONKrQiIhI5vKrLw1VaNTlFAxKaEREJGOFZjBp2nbgKaEREZHMFepyqnuWkyo0waExNCIikrlq\nBgU3NG07jveh6d+/f9zazjRKaEREJGMle9r2rbfeGre2M426nEREJHMledq2xI4SGhERyVihJ2kn\nadr2JZdcwiWXXBK39jOJupxERCRzRfMspzgmNF9//XXc2s40qtCIiEjmCnU5JWdQsMSOEhoREclY\nUQ0KdtWwpypxQcl+UUIjIiKZq5FBwaHKje5Fk/I0hkZERDJW4xWaHG+7qgpc89yYH3/AgAExbzNT\nKaEREZHM1dijD2puuBenqduTJk2KS7uZSF1OIiKSsay6wutuMqt7A7/LyarKExiV7A8lNCIikjF2\nlFexfVclzjlvQfXu+rubCKvQxGmm04UXXsiFF14Yl7YzjbqcREQkI/zpvQ3ctfDfAFzavwu/GNQd\nqnfXO2UbiHuFZteuXXFpNxOpQiMiIhlh3sebObwgl75d2vDKJ5vZ45w3KLiBCg1xHkMjsaOERkRE\n0l7xN+Ws2lzKuX0OYtQJndhaWslHG3b6FZqGupy8WU66uV7qU5eTiIikvddXfwXAoCMKyG+VTXYz\n4/XVX3FSVXloanadarqcVKFJeUpoREQk7b2+6iuOOCiXwnZe8tK/Wz6vr/qKKR3L9w78rUO8KzRD\nhw6NS7uZSAmNiIiktW1llSwr3sGEk7uGlp3e80B+/epayvPLaNVQhSY7voOCr7322ri0m4k0hkZE\nRNLaqi9LccAJh+SHlvXzfy8v39XgLKe9jz5Ql1OqU0IjIiJpbc2WUgAOL9j76IKu7VrRPMuortjV\n8BiaUIUmPgnN8OHDGT58eFzazjRKaEREJK2t2VpGu9zmHJi3dzZTdjOj+4GtcFW79o6TqUNonSo0\nKU8JjYiIpLU1W8oo6rDvgyWLCvK8wb4NDAquuUdNvCo0EjtKaEREJG3tcY61W8soKqgjoemQS3O3\nm93WwI31mmXjmmWrQhMASmhERCRtbdxeQdnuaoo65O2z7vAOubSkkm92ZzXcSFZLPZwyADRtW0RE\n0taaLWUAdXc5dcgjh938qyKLNg204bJbxu0+NCNGjIhLu5lICY2IiKStNVv3neFUo3ObLJpbNZvL\njUMbaiSrZdzuFHzllVfGpd1MpC4nERFJW2u2lNGxbUva5Oz793uz6t0AbCqzBttw2Tlxq9CUlZVR\nVlYWl7YzjSo0IiKSttZ9tYvuB7aqc13NuJiNjSQ0ZMevQnPRRRcBMH/+/Li0n0lUoRERkbRV/E05\nXQ+o5z4zfpLyVUUzyiur623DZcVvDI3EjhIaERFJSzvKq9hRXhV6IGWkmgpNuWtOyfYGEpZszXIK\nAiU0IiKSljZ84yUhXfLrqdDUJDS0CG1bF5fVUvehCQAlNCIikpZCCU09XU41d/+toDnFDSQ0XoVG\nCU2q06BgERFJSxu2e0lKYX1jaPwKjcvKCW1bl3hWaMaMGROXdjOREhoREUlLxdvKaZuTXeeUbSA0\nc6lN67xGKjQ5cavQjB07Ni7tZiJ1OYmISFoq/qa8/uoMhCo0B7Rp3fAYmuyWoW1j7auvvuKrr76K\nS9uZRhUaERFJSxu2V3DEQfveIbhGzcyldm3bsGFTBc45zOq4J00c7xR86aWXAroPTSyoQiMiImln\nj9q8a/MAABKtSURBVHOUbC+vd0AwELq3TIf8NpTuruabXVV1bxfHZzlJ7CihERGRtLNl524qq12D\nXU41FZoO7fMB6u12cjUVGudiH6jEjBIaERFJO8WN3YMGQjOXDm7XttY++8huWWt7SU1KaEREJO00\ndg8a2Fuh6dRYhSbbb0PdTilNg4JFRCTtlOzwko9ObVvWv1FVOc6yyG2VwwGtstm4o56EJctrw6or\niHWn07hx42LcYuZSQiMiImln4/YKCvKa0yK7/o4Iq6oIdSd1bNuSjfU8z2lvhSb2U7fPP//8mLeZ\nqdTlJCIiaWfTjgo6NTR+BrwKjZ+sdMpvWX+Fxk964nFzveLiYoqLi2PebiZSQiMiImln446Khrub\n8O8U7HcndWrrJTSujplMLit+g4InTJjAhAkTYt5uJlJCIyIiacU551doGk5oaldocijbXc3O8up9\nt4tjhUZiRwmNiIikla/LKqmo2kPHxio0VeVQk9D4227cse84mdAYGk3bTmlKaEREJK1simaGE0BV\nhfecJgglP3WOo8lShSYIlNCIiEhaqZmt1FiXU3iFprO/bV0znWqSnng9oFJiQ9O2RUQkrWyMtkJT\nXYFr6d1Ur11uc1pkWai6U0vYfWhi7Zprrol5m5lKCY2IiKSVjdsryG2RRduchr/irKocl3cwAM3M\n6Ni2ZeiGfOH2Vmhin9AMGzYs5m1mKnU5iYhIWtm4o4KObf9/e/cWG8d133H8+99d7pJcXiSSlihK\nrmWxUmTLlt0aiFrDhgUkQdsYaYv2oY6TIEEbpE3gB6OoixppjVZBa6B+qB/UukYltLErpYkSKS+K\nYrhAadRGreYqWaody7rYEi8yRYvk8rZcck4fZpf3yy451M5wfx9AgDk7e3zg8Wh+POd/zqQws6VP\nnLHKCfwpqp6FNtdbwxGaCxcucOHChcDbrUQaoRERkXWlZzBL23LTTczeKRhgS0M1b1y6Oe+8tXyX\n05NPPgnAyZMnA2+70miERkRE1pXugTFal9uDBuaP0DSkuDE0Tm7Sm31eYu1GaCQ4CjQiIrJujOYm\nuTk6sXxBMLN3CgZobUzhYH5hcKwKh2mVU8gp0IiIyLoxtWR7uUDj3IIjNABdc+tozCCR0j40IadA\nIyIi60axe9Dg5TDnTe1DAzP2ollopVM8pZ2CQ05FwSIism505V9d0Lbsm7b9cOJmFAW3NqQw/Bqc\neRIpfyO+gD311FOBt1mpFGhERGTd6B7IkogZt9UllzxvKpzMGKGpisdoqUvOn3ICXFUaxocD7SvA\n/v37A2+zUmnKSURE1o3ugSyb6pPEY8vsQZOfPnLx2VNTbY2phd/nlExjueADzdmzZzl79mzg7VYi\njdCIiMi60TWYnaqFWcpCIzTgFwaf6x6ad75L1mHj84+v1tNPPw1oH5ogaIRGRETWje6BMbY0LFM/\nA1NLsGfW0ABsaaymZzCL59ys464qja3BlJMER4FGRETWhdykR+/Q+PIrnGB6CfacEZq2xhQTnqN3\naHzWcZesgzUYoZHgKNCIiMi6cD0zjueKWLIN0yM0c2poCnvRdM8pDHbJujWpoZHgKNCIiMi6UFhu\n3VbElFMhnLhketbxQhiaG2hYoxoaCY6KgkVEZF0oelM9wLIZ/x9SDbOOF/avKexnU+CSab+Q2JuA\nWHCPzmeeeSawtiqdAo2IiKwLXfnl1q3FvMdp3A80Llk363htMk5jdWLBKSfAr6Op3hBAb3379u0L\nrK1KpyknERFZF7oHsrSkq0glln+0FaaPXLJ+3mdbGlPzA01VOv+9YOtoTp8+zenTpwNts1JphEZE\nRNaFzoExtm4oYsk2QDaDiychMX80Z+uGai7dGJl9MFkINEO4ed9YuQMHDgDahyYIGqEREZF14drN\nMbYVGWhsPLPg6AzA7Ruq6ewfm7UXTWHKSSudwkuBRkREIi836dEzmC0h0AzNq58p2LqhmvFJR29m\nei+a6SknrXQKKwUaERGJvM7+LA64fWORgSY7CKlFRmjybVztn17p5Arnarfg0FKgERGRyOvMh49i\na2j8EZqFA01hlOfazECjEZrQU1GwiIhE3tX+UQBu31hT1Pk2PoRXv3XBz1obUsTNr8mZskY1NM8+\n+2yg7VUyBRoREYm8qzfHqE7EaElXFXW+ZTO4loVHaKriMVobU1y9OTp1zCXXZoRm7969gbZXyTTl\nJCIikdfZ7y/ZNrPivrDEKieA2zfUzJpyIp7CxRKBv6Cyo6ODjo6OQNusVBqhERGRyLvaP1Z0QTDO\n+SM0ixQFA2zbWM1/vnNj+oAZriod+MZ6zz33HAD79+8PtN1KpBEaERGJNOcc1/qL34OGiTHMTS66\nbBv8wuD+0QkyYxPTB/WCylBToBERkUjrG84xlvNK2lQPgCWnnBZY6ZRMa2O9EFOgERGRSCvsF1N0\noMm/aXupKaetCwaa+sBraCQ4CjQiIhJpl/v89y5tby5+yTbMf9P2THc01cxq2z8/+BoaCY6KgkVE\nJNIu3xgllYjR1ljsiykHgYXftF1Qm4zT1pji8o3ppdtUpbFM92q6Os/zzz8faHuVTIFGREQi7dKN\nEbY31xCPFbdke2qEZokpJ4AdzbVcnDVCUxd4Dc3OnTsDba+SacpJREQi7VLfCO3NtUWfPz3ltEyg\naanlSt8ok57Ln58OfJXTqVOnOHXqVKBtVioFGhERiayR8Um6BrLsaCkh0OSnnBZ7OWXBjpYashMe\nXQN+YbBL1vlFwc6tuL9zHTx4kIMHDwbWXiVToBERkcgqFO3e2VJcQTAUVxQM/pQTwMUb+WmnqjTm\nPJgYW+JbUi4KNCIiElmX8mGjvZQRmvEMLlEDsaXLSO/Mt3kpXxhcqLnR5nrhpEAjIiKRdfHGKImY\nFb9LMCz72oOChuoEt9Ulp0KTq8q/oFKb64WSAo2IiETW5b4R7miqoSpewuNsfGjZguCC9pZaLuWn\ntaamqDRCE0pati0iIpF1sXeE3a3pkr5j45ll62cKdrTUcOLMdTznpqecCkXFAXjxxRcDa6vSaYRG\nREQiaWA0x9X+Me7aXFw4KSh2yglg9+Y6RnMeV/pGcelNAMSGrpfc18Vs27aNbdu2BdZeJVOgERGR\nSDrf7U/93NNWXDgpsPHMki+mnKnQ9rmuDF5dq//9oeB2Cz5+/DjHjx8PrL1KpkAjIiKR9FZXBgPu\nbi11hGaw6BGa7U01pJNx3urKQLIOl6wnNtSzgt4u7PDhwxw+fDiw9iqZAo2IiETSue4MdzbXUF9d\nQjnoZA4b/hCvbktRp8djxp4tdZzLjwZ5da2BBhoJjgKNiIhEjnOOc11DpU83DfVgzsNr2Fr0d+5p\nq+fdD4fJTnh4da2YAk0oKdCIiEjkdA9m+Wgkxz1bSgs0sUwnAF59CYFmSx0TnuMX14dw9VuIBfzG\nbQmGAo2IiETOW10ZAO5pK61+JjZwFQDXeHvR3ymMAr2VLwy2kV7wJkr698ra0z40IiISOT+7Okgq\nEWPXptL2oJkaoSmyhgZgc32KzfVJfnZtELdjM+Y8bLgXV198G4t56aWXVt2G+DRCIyIikfPGpZt8\n/I7G0nYIBmKD1/DSmyCRKul7D+7YyJuX+xlP+0u3gyoMbm5uprm5OZC2Kp0CjYiIRMqVvlE+uDnG\nQ+0bS/5ubPAaXkPpG9k91L6RTHaSd0byuwUHtBfNkSNHOHLkSCBtVToFGhERiZT/vvgRAA+3N5X8\nXRvsXFGg+fXtG0jEjNd6/JGdoEZojh49ytGjRwNpq9Ip0IiISKS8fvEmO1pq2VrCG7YBcB6xTFdJ\nK5wK0qkED/xSA69c8XDxpFY6hZACjYiIRMZQdoIffzDAwyuYbrLhDzEvh1vBCA3AQ+1NXOwbJVez\nCRsO7n1OEgwFGhERiYyT53qZ8Byf2t1S8ndjg9cAVjTlBPCJXc0Y0OOaNEITQgo0IiISCc45vvWT\nLvZsqePeEncIhpmBpvQpJ4CtG6p5ZGcT54frsYEPwLkVtSNrQ4FGREQi4c0r/VzuG+WzD7St6PuF\nTfVWUkNT8NkH2ngjt4v4UDex/ssrbqfg2LFjHDt2bNXtiAKNiIhEgHOOf3uzk6baKn7jrtKnmwAS\n77/GZPPHIFnaZnwz7dveyOWGjwMQu9yx4nYKamtrqa2tXXU7okAjIiIRcPJ8L29e6efLD24jmSj9\n0WVD14l3/ojcrkdX1Q8z43cf2cdFbwu9Z15ZVVsAhw4d4tChQ6tuRxRoREQk5HqHxvn7Vy9x/9Z6\nHlvhdFPVe6cwHLmdn151fz61u5lLDfvYfPMnXO75aFVtnThxghMnTqy6T6JAIyIiIdYzmOWPjrxF\ndtLjrx/dSTxmK2qn6t2TTDbtxGveteo+mRl7HvoMNTbO4W9/l4s3RlbdpqyeAo2IiISO5xw/OP8h\nn//mGfqGx3nhD/ZwZ/PKak3i104T7/zfVU83zZTe+TATVfU84b3MEy+/yXd+2s2Ep1VP5aS3bYuI\nSFk455jwHGM5j0x2gr7hHO9/NMr57gwdFz6iayDLxzanOfDo3ezeXFda45M5bKiHqisdVL92AG/j\nDsb3fiG4zieqyX7mBXad+BL/kniOZ1/9bY6+vpsHdm5j79YGtjfV0FyXpCGVoCYZIxEzzFY2uiTF\nMRfRdfQDAwOR63g6nWZ4eLjc3ZAi6FpFx3q8Vm8f/gr3DnSUuxtrppi/vOMxI5H/U1KrzgPnYZPj\nU59MtN7PyO/8K642+LdaV73zfWpe/XNsYgyASWdMkMDD8DAc8/s/88hvfvMmAD/8Yuk7H6+192r2\nsv2r3y13N2hsbCzqf4LIBhoRERGRAtXQiIiISOQp0IiIiEjkKdCIiIhI5CnQiIiISOQp0IiIiEjk\nKdCIiIhI5CnQiIiISOQp0NxCZrbTzMbM7N/nHH/czN43s2Ez+76ZNZWrj5XMzFJmdjh/LTJm9nMz\n+60553zCzN4xsxEz+y8zu6Nc/a10ZtZkZify9837ZvZ4ufskvuXuJd1H4bTQMypKzycFmlvrH4Ef\nzTxgZnuAF4EvAJuBEeCfbn3XBP9VIFeBR4BG4C+B75jZdgAzawGOA38FNAE/Br5djo4K4N9P4/j3\nzeeAF/L3k5TfoveS7qNQm/WMitrzSTsF3yJm9hjwe8D/Ab/snPt8/vjfAdudc4/nf24H3gaanXOZ\ncvVXfGZ2Fvgb59z3zOwrwJeccw/mP0sDN4Bfcc69U85+Vpr8f/ubwD3OuXfzx14GOp1zf1HWzsmC\nCvcS0Izuo9BZ6BkVteeTRmhuATNrAA4Af7rAx3uAM4UfnHMX8X/rXP077mVVzGwz/nU4nz8091oN\nAxfzx+XW2gVMFMJM3hl0LUJpzr2k+yhklnhGRer5pEBza3wDOOycu7bAZ3XAwJxjA0D9mvdKFmVm\nVcAR4JszfmvUtQqPOmBwzjFdixBa4F7SfRQ+iz2jInWtFGhWycw6zMwt8ud1M7sf+CTwD4s0MQQ0\nzDnWAIRuOC/qlrtWM86LAS/j/ybyxIwmdK3CQ9ciAha5l3TtQmSZZ1SkrlWi3B2IOufc/qU+N7Mn\nge3AB2YGfuKNm9ndzrlfxR+CvW/G+TuAFPDu/NZkNZa7VgDmX6TD+AVwn3bO5WZ8fB744oxz00A7\n01NScuu8CyTMbKdz7kL+2H3oWoTGEveS7qNw2c8izyjgh0To+aSi4DVmZrXMTrh/hv8/z1edc735\nKvL/AR4FfopfUZ5wzj12q/sqYGb/DNwPfNI5NzTns9uA94A/BE7iFzg+4pz7tVveUcHM/gNwwJfx\nr9kPgAedc3owhsBi95Luo3BZ6hkFbCJCzydNOa0x59yIc66n8Ad/CG/MOdeb//w88Cf4c8wf4s9N\nfq1sHa5g+b0w/hj/L+EeMxvK//kcQP6a/T7wt/grbPYBobyxK8TXgBr8++Zb+L8kKMyEwFL3ku6j\ncFnqGRW155NGaERERCTyNEIjIiIikadAIyIiIpGnQCMiIiKRp0AjIiIikadAIyIiIpGnQCMiIiKR\np0AjIiIikadAIyIiIpGnQCMioWZmX89vo1/4eaOZ5cysupz9EpFwUaARkbC7F/j5jJ/vB37hnBsr\nU39EJIQUaEQk7BYKNGfK1BcRCSkFGhEJLTNLAu3A2RmH72N2wBERUaARkVC7C+h0zo0AmJkB+9EI\njYjMoUAjImG2F9hkZu1mVgN8A7gDuFLWXolI6CjQiEiY3Qu8AnQA7wEZ4Brw9TL2SURCyJxz5e6D\niMiCzOwUcMg5971y90VEwk0jNCISZvcCb5e7EyISfhqhEZFQMrONwHUg7ZzLlbs/IhJuCjQiIiIS\neZpyEhERkchToBEREZHIU6ARERGRyFOgERERkchToBEREZHIU6ARERGRyFOgERERkchToBEREZHI\n+38btR4mLesAdgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2a3098c160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plotting\n",
    "\n",
    "# create figure\n",
    "fig, axes = plt.subplots(2, 1, sharex=True, figsize=(8, 8))\n",
    "\n",
    "# plot histogram\n",
    "ax = axes[0]\n",
    "ax.hist(y, np.arange(-44, 43, 2))\n",
    "# decorate\n",
    "ax.set_title('Newcomb\\'s measurements')\n",
    "ax.set_ylabel('count')\n",
    "ax.set_xlabel('$\\mu$')\n",
    "plt.setp(axes[0].get_xticklabels(), visible=True)\n",
    "\n",
    "# plot the posterior of mu\n",
    "ax = axes[1]\n",
    "ax.plot(t1, pm_mu)\n",
    "# plot the posterior of mu in the filtered case\n",
    "ax.plot(t1, pm_mu_pos)\n",
    "# Plot the currently accepted true value\n",
    "ax.axvline(33, color='k', linestyle='--')\n",
    "ax.legend(\n",
    "    ('posterior of $\\mu$',\n",
    "     'posterior of $\\mu$ given $y > 0$',\n",
    "     '\"true value\"'),\n",
    "    loc='upper left'\n",
    ")\n",
    "ax.set_title('Normal model')\n",
    "ax.set_xlabel('$\\mu$')\n",
    "ax.set_yticks(())\n",
    "# set bottom to zero\n",
    "ax.set_ylim((0, ax.set_ylim()[1]))\n",
    "\n",
    "fig.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
